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The steady, coaxial flow in which two immiscible, incompressible fluids move past each 
other in a cylindrical tube has a continuum of possibilities due to the arbitrariness of the 
interface between the fluids. By invoking the presence of surface tension to at least restrict 
the shape of any interface to that of a circular arc or full circle, we consider the following 
question: which flow will maximise the exchange when there is only one dividing interface 
r? Surprisingly, the answer differs fundamentally from the better-known co-directional 
two-phase flow situation where an axisymmetric (concentric) core-annular solution always 
optimises the flux. Instead, the maximal flux state is invariably asymmetric either being 
a 'side-by-side' configuration where T starts and finishes at the tube wall or an eccentric 
core-annular flow where T is an off-centre full circle in which the more viscous fluid 
is surrounded by the less viscous fluid. The side-by-side solution is the most efficient 
exchanger for a small viscosity ratio f3 < 4.60 with an eccentric core-annular solution 
optimal otherwise. At large /3, this eccentric solution provides 51% more flux than the 
axisymmetric core-annular flow which is always a local minimiser of the flux. 



1. Introduction 

For Newtonian fluids at least where the governing Navier-Stokes equations are known, 
the most fundamental issue in fluid mechanics is predicting the realised flow solution for 
a given initial state and set of boundary conditions against a background of omnipresent 
noise. Non-uniqueness of solution is endemic due to the nonlincarity of the Navier-Stokes 
equations but even in special limits (e.g. vanishing Reynolds number or steady, unidi- 
rectional flow) where these simplify to the linear Stokes' equations, degeneracy is rife as 
specification of the flow domain is typically part of the problem. A well-known exam- 
ple of this is the pressure-driven flow of two immiscible fluids along a cylindrical tube 
(e.g. Joseph, Renardy & Renardy 1984, Joseph, Nguyen and Beavers 1984, and Joseph 
et al. 1997). Here there is a continuum of steady unidirectional solutions possible due 
to the arbitrariness in the interface between the two fluids. In practice, however, the 
axisymmetric core-annular solution with the more viscous fluid surrounded by the less 
viscous fluid is invariably observed for fluid combinations ranging from oil and water 
(Charles & Redberger 1962, Yu & Sparrow 1967, Hasson, Mann & Mr 1970), to molten 
polymers (Southern & Ballman 1973, Everage 1973, Lee & White 1974, Williams 1975 
and Minagawa & White 1975). 

Interestingly, it appears that if an extra constraint is added to the system - that the 
mean volumetric flux along the tube vanishes - different steady solutions are observed 
(Arakeri et al. 2000, Huppert & Hallworth 2007, Beckett et al. 2009). Such a flow is 
easily set up in the laboratory by placing a tank of dense fluid directly above a tank full 
of less dense fluid and connecting the two by a vertical cylindrical tube. If the density 
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difference or the tube cross-section is small enough or the fluid viscosities large enough, 
it is reasonable to anticipate a steady, coaxial flow established in the tube in which the 
denser fluid falls under gravity displacing the less dense fluid upwards. When the lower 
tank is initially full and both fluids incompressible, this exchange flow is constrained to 
have no net volume flux along the tube. As in the unidirectional flow situation, the form of 
the steady, coaxial two-fluid flow realised is fascinatingly unclear due to the arbitrariness 
of the interface between the fluids (formally, any union of open curves terminating on 
the tube wall and closed curves in the interior are possible). Using salty and pure water, 
Arakeri ct al (2000) saw only a 'half-and-half solution where the interface divides the 
tube cross-section into two approximately equal domains (hereafter referred to as a 'side- 
by-side' solution). In contrast, Huppert & Hallworth (2007) saw only a concentric core- 
annular flow as their steady low-Reynolds solution and recently both types of flow have 
been seen in the same apparatus (Beckett et al. 2009). Beyond its intrinsic interest, 
this flow has applications ranging from the exchange of degassed and gas-rich magma in 
volcanoes (e.g. see Huppert & Hallworth 2007 and references herein) to plug-cementing 
oilfields (e.g. Frigaard & Scherzer 1998, Moyers- Gonzalez & Frigaard 2004). There is also 
associated work on exchange problems involving misciblc fluids, tilted tubes or channels, 
and unsteady solutions (see the recent articles by Seon et al. 2007, Znaien et al. 2009 
and Taghavi et al. 2009 for references). 

Resolving the flow degeneracy of the steady state in favour of one realised solution 
involves knowledge of the initial conditions of the exchange flow, the pressure boundary 
conditions set-up across the tube and the inherent instability mechanisms present. Prag- 
matically, the initial conditions are never known that well (e.g. barriers are slid open 
or plugs removed in the laboratory), the pressure gradient which gets set up difficult 
to measure and assessing relative stability requires every possible flow state to be iden- 
tified first. It is therefore tempting to jump to an ad-hoc selection principle especially 
as a particularly obvious one suggests itself here: the flow selects the solution which has 
the largest individual volumetric flux. A selection principle based upon maximum flux 
has some history in the undirectional two-phase flow problem motivated by its formal 
connection to the single fluid problem (Maclean 1973, Everage 1973, Joseph, Nguyen & 
Beavers 1984). Here, the governing Stokes equations are the Euler-Lagrange equations 
for maximising the flux for velocity fields which satisfy the global power balance that 
the rate at which energy is viscously dissipated equals the power supplied by the applied 
pressure gradient (per unit length of the tube). Specifically, if G is the constant applied 
pressure gradient, f2 the cross-section of the tube and u the speed along the tube, then 



where 8 indicates the Frechet (variational) derivative, J —GudA is the rate of working 
by the pressure gradient per unit length of tube and the Lagrange multiplier A imposing 
the power balance constraint takes the value 1/G. The stationary point defined by the 
variational solution is clearly one of maximum flux because the only quadratic term in the 
integrand is negative definite (u is oppositely signed to G so A < Ojjj. The fact that this 
variational formulation can be extended to two fluids provided the interface between them 
is known (Maclean 1973, Everage 1973) supplied the impetus to invoke the principle of 
maximal flux more generally. It appears to be mostly successful - in the words of Joseph, 

| Due to the relative simplicity of Stokes equations, there are many other variational formu- 
lations such as maximising the dissipation subject to the global power balance, minimising the 
dissipation subject to fixed flux and the complementary problem of maximising the flux subject 
to fixed dissipation. 





(1.1) 
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Nguyen and Beavers (1984) "our experiments show that something like this is going on"- 
predicting that the more viscous fluid will be encircled by the less viscous fluid which 
then acts as a lubricant against the tube walls (see also Charles & Redberger 1962, Yu & 
Sparrow 1967, Hasson, Mann & Nir 1970, Southern & Ballman 1973, Everage 1973, Lee 
& White 1974, Williams 1975, Minagawa & White 1975). Joseph, Renardy & Renardy 
(1984), however, add some qualifications: this state can become unstable if the more 
viscous core gets too small. 

Given this history, the purpose of this paper is to explore the consequences of this 
'maximum flux principle' in predicting the form of the exchange flow realised in a ver- 
tical cylindrical tube. Formally solving the variational problem with the interface (or 
interfaces) as an unknown is a formidable challenge not attempted here. Rather, a sur- 
vey is conducted over a physically-motivated subspace of all mathematically-possible 
steady, coaxial solutions. This subspace is defined by two (mild) assumptions: a) the 
fluids occupy one (possibly multi-connected) domain so that there is only one interface 
r, and b) that this interface is a circular arc or a full circle. The motivation for the for- 
mer assumption is stability - multiple small fluid domains would presumably aggregate - 
and the presence of some surface tension between the two fluids conveniently motivates 
the latter. The axially-constant, lateral pressure difference required to balance interfa- 
cial tension, however, will be ignored in what follows as it has no consequence for the 
calculations. 



2. Formulation 

Consider two immiscible fluids with densities p\ and p 2 and viscosities p\ and p 2 which 
are flowing in a vertical circular tube of radius a across which there is a pressure gradient 
G and g is the acceleration due to gravity. Assuming that fluid 1(2) occupies an area 
A\{A 2 ), the Navier-Stokes equations for steady exchange flow of the two fluids either 
directed up or down the tube (so the problem is just in the cross-sectional plane) are 

G = p{W 2 u\ - pig in A\, G = p 2 \7 2 u 2 - pig in A* 2 (2.1) 

with non-slip boundary conditions at the tube wall and continuity of velocity and stress 
at the interface T* between the two fluids, that is 

u * =u * & Ml M = Al2 £g on r * (2.2) 

(where d/dn is the normal derivative to T*). There is a further constraint that the net 
volume flux through the tube is zero so 

Q*:=-J u\dAl= j u* 2 dA* 2 . (2.3) 

Without loss of generality, we assume p\ > p 2 so that Q* is positive (the less dense fluid 
rises). This does not prejudice the choice of viscosities later because of the symmetry 
(Pi>P2>5) —* (P2,Pi,—g)- the direction 'up' is irrelevant with only the density difference 
being important. 

The system is non-dimensionalised (*'s removed) using the tube radius a, the differ- 
ential hydrostatic pressure gradient Apg (where Ap := pi — p 2 ) and pi so that after 
defining A by 

G = -±( Pl + p 2 )g + ±Apg\ (2.4) 
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Figure 1. The side-by-side solution configuration specified by two parameters: 7 and a. 
then 

VV=A + 1 in A t , (2.5) 
/3V 2 u 2 = A-l in A 2 , (2.6) 

dut du 2 

ui = u 2 & — — = p— — on T. (2.7) 
on on 



where 

El 
mi ' 

Henceforth u\ and u 2 are in units of ^Apga 2 / [i\ and the one-fluid volume flux 



/?:= 



Q'=-J u 1 dA 1 = J u 2 dA 2 (2.9) 

is in units of ^Apga 4 /ui with Ai U A 2 being the unit disk. 

Two specific choices are now made for T. The first is a circular arc of general curvature 
and position which intersects the tube wall so that the two fluids are next to each other 
- the side-by-side solution: see figure [TJ The second is a full circle completely contained 
within, but not concentric with, the tube so that one fluid encapsulates the other - the 
eccentric core-annular solution: see figurc[2l The limiting case of a concentric core-annular 
solution needs to be treated separately but is easily solved analytically. 



2.1. Side-by-side solutions 

The geometry of the side-by-side solution is shown in figure Q] to be defined by two 
parameters: 7, the (upper) intercept latitude of T with the tube wall, and 2a, the angle 
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Figure 2. The eccentric core-annular configuration specified by two parameters: a and R. 

between T and tube wall. For given viscosity ratio (3 and pressure gradient A, one of 
these (nominally a) is determined by the flux balance leaving a 1-dimensional family of 
side- by-side flows with corresponding fluxes Q = Q s (f3, A; 7) possible (see appendix A for 
the calculation details). There is a symmetry 

Q(/3, A; 7, a) = 4<3(4, -A;tt - 7, | - a) (2.10) 

which means that only [3 ^ 1 need be considered providing the full ranges of 7 and a 
are studied. Henceforth fluid 2 will always be the more viscous fluid so that the non- 
dimensionalisation has been done using the smaller dynamic viscosity Hi. 



2.2. Eccentric solutions 

The eccentric core-annular solution has one fluid domain as a totally-contained circular 
disk (cylinder) not touching the tube wall. The radius R < 1 and centre (a, 0) of T define 
the geometry uniquely up to obvious rotations and reflections. To match smoothly onto 
the choices made in the side-by-side solution, a is chosen to be +ve(-ve) for A\ in A 2 
(A 2 in Ai). As before, for given viscosity ratio (3 and pressure gradient A, one of these 
two geometrical parameters is determined by the flux balance. This is done by searching 
over R for given 

, _ / 1 + cr - R 4 2 in 4i cr < , , 

\ -l + a + R Ai in 4 2 a > { ' 

which either represents the positive displacement from (—1,0) to (cr — R, 0), the leftmost 
point of r for the case of A 2 in A\ (a < 0) , or the negative displacement of (a + R, 0), 
the rightmost point of F, from (1, 0) for the case of A\ in A 2 (a > 0). This choice is made 
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for two reasons. Firstly, d is a convenient way of extending the side-by-side solutions 
continuously beyond their pinch-off points into the corresponding eccentric solutions: 
7^0 corresponds to A2 encapsulating Ay and d decreasing across zero whereas 7 — > 7r 
corresponds to A\ encapsulating A2 and d increasing across zero (see figure 3). Secondly, 
only one flux-balanced solution was ever found for a given d whereas some a can have two 
flux-balanced solutions. The result is that two 1-dimensional families of eccentric core- 
annular flows with corresponding fluxes Q e (f3,A;d) (more viscous core) and Q e (ft,\,d) 
(less viscous core) are possible (see appendix B for the calculation details). It's worth 
re-emphasizing here that ft ^ 1 so all the flux values quoted are in units of where 
jUi is the smaller dynamic viscosity. 

2.3. Concentric solutions 

When r is a circle concentric with the tube wall there is a simple solution to the problem 
d23])-(!27]) discussed recently by Huppert & Hallworth (2007): 



«2 



^±V-l)-i? 2 logr, 

A- 1 

4/3 



(r 2 -R 2 )-R 2 logi? 



A+l 



R < r 1 



(1-i? 2 ). rs^R 



The associated fluxes are 



^- 8/3 



(A + l)(2i? 2 - i? 4 - 1) + 4i? 2 (l - R 2 ) + 8i? 4 logi? 
(1 - A)i? 4 - 2/3(1 + A)i? 2 (l - R 2 ) - 8ftR A logi? 



Since this is a special case of an eccentric core-annular solution with a 
unique < R < 1 for a flux-balanced solution which is 



'2/3 - V4/3 2 - (3(1 + A)[/3(3 - A) + (A - 1)] 



/ [/3(3 - A) + (A - 1)] 

so that the flux (for fluid 2 in the core) is Q c (ft, A). As f3 — » oc 
R c -> y/(l + X)/(3-\), A -> 0.1746 and Q c 



0.01831 



(2.12) 
(2.13) 

(2.14) 

(2.15) 
0, there is 

(2.16) 
(2.17) 



from above. The opposite scenario of the less viscous fluid (fluid 1) in the core has 
Q := Q c - 0(/3) (/3 -> in expressions ([2~T4"1) and (|2T5|) and multiply Q by 1/0 to 
convert the flux units to those using the smaller dynamic viscosity). 

2.4. Strategy 

The strategy now is to calculate max\Q as a function of ft over all possible geometries 
smoothly ranging from the concentric solution with less viscous fluid in the core through 
to the concentric solution with the more viscous fluid in the core. Figure 3 illustrates 
the spectrum of possibilities and a glimpse of how the flux varies at one ft value. Before 
detailing the results further, the reader may be amused by an admission. At onset, this 
author (naively?) expected the calculation of maximum flux to be a simple competition 
between a local maximum achieved by the side-by-side solution and the flux Q c associated 
with the concentric core-annular flow influenced by the known behaviour of unidirectional 
2-fluid flow. The side-by-side solution, however, quickly loses its interior maximum (0 < 
7 < tt) as ft increases in favour of an end-point maximum at 7 = n. The fact that this 
end-point maximum exceeds the concentric solution flux Q c unequivocally indicated the 
importance of the intermediate eccentric core-annular flux Q e . 
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Figure 3. The flux maxAQ plotted across the various flow configurations (d < indicates less 
viscous core and d > more viscous core) for (3 = 5. The leftmost point is Q c , beyond this, the 
region d < is the domain for Q e , the region 7 6 [0, n] is the domain for Q s , d > the domain 
for Q e and the rightmost point is Q c . The interior local maxima are highlighted with dots. The 
curve is only C° because the abscissa changes character at 7 = and n of course. 



3. Results 

There is a special case of the problem which can be solved using known results. When 
(3 = 1, the optimal balanced flow of fluid 1 must mirror that in fluid 2. In particular, 
A = 0, r is the diameter x = and u = on T. The problems for either fluid then decouple 
into single phase pressure-driven flow in a 'half-cylinder (semicircular cross-section). The 
flux is 0.07438920 in our non-dimensional units according to White's (1991) equation (3- 
44) . This provides an excellent test of the side-by-side computations (see Table 1 which 
shows 3 significant figure correspondence although there is really 5). Further checks are 
available between the very different side-by-side and eccentric flow codes (e.g. figure 4 
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Figure 4. Plotting max^Q against d at j3 = 5 for side-by-side solutions as d — > 0~ (7 — > 7r) and 
eccentric solutions as d — > + demonstrates the smooth connection between the two formulations. 
Velocity fields for the circled points are shown in figure [5] 




Figure 5. The pinching-off side-by-side (d = —0.0198, —0.051 < u < 0.041) and near-touching 
eccentric solutions (d = 0.0205, -0.052 < u < 0.043) for f3 = 5 and the optimal A = -0.20 
corresponding to the circles in figure [4] The contours range from -0.105 (dark/red) to 0.105 
(light/white) in steps of 0.01 here and throughout figures [7] and [8] to aid comparison. 
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Figure 6. max^Q as a function of 7 for the side- by-side solutions and as a function of d > 
for the eccentric solutions at /3 — 2, 4, 5 and 10. The single dot at the right end of each curve 
corresponds to the concentric case Q c - The global flux maximum is a side- by-side solution for 
j3 ^ 4.60 and an eccentric solution for (3 ^ 4.60. 



where using d as the abscissa shows at least C 1 continuity in max a Qat7 = 7ror<i = 
at [3 = 5). 

3.1. max\Q 

The /3 value chosen in figure [3] has been purposely chosen to show the presence of flux 
maxima in the side-by-side solutions and the (d > or more viscous fluid in the core) 
eccentric solutions (Q e )- The complementary eccentric solutions with the less viscous 
fluid in the core (Q e ) always show monotonic behaviour in which the flux decreases from 
the 7 = side-by-side value down to the concentric core-annular value of Q c (leftmost 
point or most negative d). This uninteresting part of the flux spectrum is suppressed in 
figure[6] to focus on max\Q over 7 and d > for (3 G [2, 10] over which all the interesting 
behaviour occurs. At f3 = 1, the side-by-side solution with 7 = n/2 and a — 7r/4 supplies 
the only flux maximum with both concentric core-annular solutions being global minima 
as Q c — Q c . At (3 « 2, a local maximum starts to appear in the eccentric solutions 
with d small and positive (see figure[5]). At /3 « 4.60, this 'eccentric' maximum becomes 
the global maximum with the 'side-by-side' local maximum disappearing by [3 w 8.2. 
Thereafter the sole flux maximum is always an eccentric solution. Figures [7] and [5] show 
how the maxima change with /3 including an eccentric optimal flux solution at (3 — 10, 000. 
This confirms that the optimal asymptotic solution has plug flow for the more viscous 
core. Figure [9] plots the maxima values as a function of (3 highlighting the cross-over 
point at j3 « 4.60 (see also Tables 1 and 2). The concentric core-annular flux values for 
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Figure 7. Maximal flux side-by-side solutions for j3 = 1 (top left, -0.098 < u < -0.098), = 2 
(top right, -0.078 < u < 0.063), /3 = 5 (bottom left, -0.058 < u sC 0.037) and (5 = 8 (bottom 
right, -0.047 ^ u ^ 0.029). The contours range from -0.105(dark/red) to 0.105 (light/white) in 
steps of 0.01 here and throughout figures [5] and [8] to aid comparison. 



the more viscous fluid in the core Q c and less viscous fluid in the core Q c are also shown 
as a local and global minima respectively. 

3.2. max\Q for f3 — > oo 

At large /3, there is every reason to suspect that the maximal flux possible possesses a 
simple expansion around its limiting value: 



maxx,dQe((3, A; d) = Qoo + + -| + . . . 

P P 



The scalars Qoo and a\ can be estimated as follows 



Qc 



Pi -I3 2 



a 1 



P2-P1 



Q(Pi) - Q(fo) 



(3.1) 



(3.2) 
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p 


A 


7 


1 


0.00 


1.57 


.5 


-0.06 


1.52 


2 


-0.10 


1.50 


.5 


-0.13 


1.48 


3 


-0.15 


1.46 


.5 


-0.18 


1.46 


4 


-0.19 


1.47 


.5 


-0.21 


1.47 


5 


-0.22 


1.48 


6 


-0.25 


1.52 


7 


-0.27 


1.58 


8 


-0.28 


1.69 



a Q s (xl0~ 2 ) 

0.785 7.44 

0.810 6.11 

0.812 5.34 

0.814 4.82 

0.813 4.43 

0.809 4.13 

0.786 3.88 

0.779 3.67 

0.761 3.50 

0.721 3.21 

0.668 2.98 

0.585 2.79 



Table 1. maxA, 7 Q s (Q for the side-by-side solution) as a function of (3. The maximum is 
unique global for (3 < 4.60 and thereafter is a local maximum until it vanishes for a (3 w 8.2. 



p 


A 


a 


R 


Qe (x: 


2.5 


-0.140 


-0.393 


0.594 


4.25 


3 


-0.155 


-0.390 


0.590 


4.02 


3.5 


-0.168 


-0.387 


0.584 


3.86 


4 


-0.180 


-0.387 


0.582 


3.74 


4.5 


-0.188 


-0.386 


0.578 


3.65 


5 


-0.198 


-0.386 


0.574 


3.58 


6 


-0.209 


-0.385 


0.570 


3.47 


7 


-0.211 


-0.382 


0.568 


3.40 


8 


-0.222 


-0.383 


0.565 


3.34 


9 


-0.225 


-0.383 


0.564 


3.29 


10 


-0.226 


-0.381 


0.563 


3.26 


15 


-0.246 


-0.383 


0.555 


3.15 


20 


-0.245 


-0.381 


0.556 


3.10 


50 


-0.260 


-0.381 


0.549 


3.01 


100 


-0.260 


-0.380 


0.551 


2.98 


200 


-0.263 


-0.381 


0.550 


2.96 


500 


-0.263 


-0.3795 


0.5485 


2.9544 


1000 


-0.263 


-0.3794 


0.5476 


2.9514 


2000 


-0.263 


-0.3797 


0.5473 


2.9499 


5000 


-0.263 


-0.3795 


0.5479 


2.9490 


10000 


-0.263 


-0.3795 


0.5479 


2.9487 


oo 


-0.263 


-0.3795 


0.5480 


2.9484 



Table 2. maxA,d<2e (Q for the eccentric solution) as a function of f3. The maximum appears 
for (3 ~ 2, is a local maximum for 2 < (3 < 4.60 and becomes a unique global maximum for 
13 > 4.60. 
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Figure 8. Maximal flux eccentric core-annular solutions for j3 = 2.5 (top left, 
-0.059 < u < 0.058 ), = 5 (top right, -0.053 < u < 0.044), /3 = 8 (bottom left, 
-0.051 < u < 0.039) and j3 = 10,000 (bottom right, -0.047 < u < 0.031). The contours 
range from -0.105(dark/red) to 0.105 (light/white) in steps of 0.01 here and throughout figures 
[S] and [7] to aid comparison. 



where (3\ and fa have suitably large values. There is good evidence that Qoo ~ 2.9484 x 
10~ 2 and ai « 3.00 x 10~ 2 supporting the original assumption: see figure ITTJl Another 
check on this value of is available by artificially imposing plug flow in the core (e.g. 
see the lower right solution in figure [5]) . The matching conditions at V then simplify 
to just continuity u% = U2 and the condition that the continuation of u\ into Ai has 
no logarithmic singularities (^ r dx.Vwi = 0) which eliminates (3 from the problem. A 
straightforward search over A and a then reveals the maximum of = 2.94844 x 10 -2 
at A = -0.263 a = -0.3795, R = 0.54798 (and d=l + a-R = 0.0725). 
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Figure 9. maxA^Qs (left upper black curve), maxA,dQe (right upper blue curve) and max^Qc 
(lowest solid red curve) compared as a function of (3. The side-by-side maximum disappears for 
/3 > 8.2 and the eccentric solution only starts to have a maximum for f3 > 2. The lowest dashed 
(red) curve corresponds to max\Q c , the global minimum of the more viscous fluid encapsulating 
the less viscous solution. 



3.3. Q ((3, X) for fixed X 

So far all the results shown have been optimised over the pressure gradient A. The 
presumption is that, in the absence of any explicitly imposed gradient, the flow sets up 
its own to maximum the volumetric exchange. Figure [IT] shows the effect of fixing A on 
the flux profile at /3 = 5. The same general trends emerge with one important additional 
feature highlighted by the A = —0.5 curve. Here Q c (leftmost point) is approximately 
the same as Q c (rightmost point) . Figure [12] plots the two core-annular flux functions 
Q c and Q c against A to show that the less-viscous core solution flux Q c actually exceeds 
the more-viscous core solution flux Q c for A < —0.51 at (3 — 5. This threshold pressure 
gradient montonically decreases as (3 increases to, for example, k, —0.89 at (3 — 100 (recall 
— 1 < A < 1): see figure IT2"1 Since a A value of -1 translates into a pressure gradient which 
hydrostatically maintains the denser fluid, the conclusion is that the less-viscous-fluid- 
in-the-core concentric solution is favoured over its complement for large enough pressure 
gradients. 



4. Discussion 

This paper has considered the steady, coaxial flow of two immiscible fluids of different 
densities and viscosities in a straight vertical cylindrical tube such that their volumetric 
fluxes balance. Under mild assumptions concerning the interface between the two fluids, 
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Figure 10. The coefficients Qoo and a\ (inset) as calculated using the expressions (|3.2p 
against j3 where j3 and the next smallest value of /3 were used. 

the main conclusion is that the flow which optimises the volumetric flux over all possible 
pressure gradients is always asymmetric. In particular, for viscosity ratios < 4.60 the 
optimal flow is a side-by-side solution in which each fluid makes contact with a side of 
the tube and otherwise is an eccentric core-annular solution with the more viscous fluid 
encapsulated by the less viscous fluid. (In fact, in this latter case, the eccentricity is so 
marked, that it could look like a side-by-side solution from one direction to the unwary.) 
The axisymmetric (concentric) core-annular solution in which one fluid encircles the other 
is surprisingly either a local or global minimiser of the flux. The clear conclusion is that 
displacing the core of such a flow to one side increases the flux by allowing the outer fluid 
to 'bulge' through the larger gap. This generalises the equivalent observation made for 
the flow of a single fluid through an eccentric annulus duct (see figure 3-8 on page 127 
of White 1991). 

The fact that the principle of maximum flux predicts a side-by-side solution at low 
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Figure 11. The effect of fixing the pressure gradient A at —0.5, and 0.5 on Q for /3 = 5. The 
(black) dashed upper envelope is the result of optimising over A as shown in figure [6] 



viscosity ratios does find support in the work of Arakeri et al (2000) and the experiments 
at Bristol (Beckett et al. 2009). However, Huppert & Hallworth (2007) never mention see- 
ing a side- by-side solution during their low- viscosity-ratio experiments, instead reporting 
only a steady concentric core-annular flow. More intriguing, however, is that in this core- 
annular solution, both Huppert & Hallworth (2007) and Beckett et al. (2009) invariably 
see the lower (less dense) fluid rising along the axis. On the basis that less dense fluids 
generically are less viscous too, this implies that the less viscous fluid is typically at the 
core of these observed flows. From the flux perspective, the results presented here show 
that this globally minimises the flux over all possible pressure gradients! This apparent 
contradiction is ameliorated somewhat if the pressure gradient set up (or imposed) is 
towards the maximum possible for exchange (e.g. see figure IT2|). but nevertheless the 
core-annular solution still remains a local flux minimiser. The principle of minimum flux 



16 



R. R. Kerswell 



0.025 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



X 

Figure 12. The concentric core-annular fluxes Q c (/3, A) and Qc(f3, A) plotted against A for /3 = 5 
(thick solid red and thin solid black respectively) and f3 = 100 (thick dashed red and thin dashed 
black respectively). The crossing of the solid lines at ~ —0.51 is consistent with figure [TT] where 
Qc « Qc at A = —0.5. The dashed lines cross at A m —0.89 for a ratio of 100. 

(and, coincidentally, minimum dissipation) then appears more useful at large viscosity 
ratios. 

The proper route to resolving this conundrum, of course, is careful consideration of the 
initial value problem and the stability of the evolving solution to the small disturbances 
always present. A first step in this direction would be to study the Rayleigh- Taylor 
instability problem in a cylindrical tube where a fluid of density p\ and viscosity pi fills 
the half cylinder z > and a fluid of density p2 < pi and viscosity pi occupies z < 0. 
Establishing which interfacial deformation mode (axisymmetric or asymmetric) has the 
largest growth rate as a function of all the parameters present would surely go some way 
in predicting which type of flow is initiated. However, even this calculation doesn't seem 
to have been done yet although Batchelor & Nitsche (1993) come close. 
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In conclusion, it should be clear that there are some interesting issues surrounding the 
exchange flow of two fluids in a vertical tube. Even the steady immiscible problem displays 
an intriguing degeneracy of solution. Focussing on an ad hoc principle of maximum (or 
minium) flux unfortunately looks to be too simplistic despite its appealing rationale and 
apparent success in an associated context. This means that there is no avoiding a more 
formal stability-based approach to explain what is seen in experiments. 
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Appendix A. Side-by-side solutions 

The geometry of the side-by-side solution is shown in figure [T] to be defined by two 
parameters: 7, the (upper) intercept latitude of T with the duct wall, and 2a, the angle 
between T and duct wall. The coupled Poisson problems (|2.5[) - (|2.7p become two Laplace 
problems by separating off simple inhomogeneous parts as follows 

ul = <f> 1 + ±±±(x 2 +y 2 -l), u * = <f> 2 + ^-(x 2 +y 2 -l) (Al) 

which have been designed to leave the boundary conditions on the duct wall undisturbed. 
The functions $1 and $2 then satisfy 

V 2 $i = in Ai, (A 2) 

V 2 $2 = in A 2 , (A3) 

with boundary conditions 

$1=0 on x + iy = e v ^ - 7 s$ -0 < 7 (A 4) 

$2 = on x + iy = e 1 ^ 7 < ip < 2tt - 7 (A 5) 

Here the interface curve T:={z\z = x + iy = <r + Re ; \0\ ^ \0 max :— 7 — 2a\ } where 

sin 2a „ „ sin 7 , . _. 

<?■■=—„ & R:= . / (A 7) 

Sill Umax Sill Umax 

are formulae for the centre (x,y) = (a, 0) and radius of curvature respectively valid for 
any pair ^ 2a, 7 $J n. (The singular case 7 = 2a where R — > 00 so that T is a straight 
line cannot be formally handled but is never a practical problem.) 

The solution strategy is to transform regions A\ and A2 into the upper and lower half 
planes respectively via conformal tranformations where an explicit solution can then be 
deduced by Poisson's integral formula. Three simple transformations prove sufficient, the 
first is 

z — cos 7 , . „. 

w := (A 8) 

sin 7 
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Figure 13. The composite conformal mapping from z to £. The interface curve V meets the 
duct wall at e ±l7 in the z-plane, which are moved to ±i in the iu-plane, and then to ±ioo in 
the q-plane. Two separate transformations then map the strip ^7 — a ^ Re(q) ^ ^7 onto the 
upper half plane and the strip 57 — tv ^ Re(q) ^ ^7 — a onto the lower half plane. In both 
cases, Im(q) > in mapped inside the appropriate unit semicircle in the £-plane. 



which rescales the duct so that its radius becomes l/sin.7 and the 2 contact points 
of r with the duct wall e ±I7 move to ±i (other noteworthy images are: — > — cot 7, 
1 — > sin7/ (1 + COS7), — 1 — » — sin7/(l — COS7) ; see Figure 2). A second transformation 

q = r + is := tan -1 w r,s £ Se (A 9) 

converts all the circular arcs into straight lines parallel to the imaginary axis in the 
complex q plane. To see this, consider the transformation in reverse 

e 2iq - 1 

w = tanq=-i e2 . q + i (A 10) 

and decompose this transformation into its 3 components. A strip ^7 — a ^ r ^ ^7 with 
< a ^ 7r is rotated through tt/2 by q — > iq. Doubling and exponentiating q — > iq — > e 2j<? 
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then tranforms the strip into the interior of a wedge centred at the origin with sides of 
argument 7 — 2a and 7. Finally the Mobius transformation q —> iq — > e 2lq — » —i{e 2lq — 
l)/(e 2lq + 1) converts the wedge sides into circular arcs joining the points w = H and 
the wedge interior into a circular lune with angle 2a (see Figure 2 and pages 205-207 of 
Marushevich 1965). The conformal transformation (|A 9p is undoubtedly not the only one 
which would do the job (e.g. Vlasov 1986) but is particularly nice since it can used to 
treat both 'lunes' together: A\ maps to the strip ^7 — a ^ r ^ ^7 and A 2 maps into the 
strip Jj7 — i-7r ^ r ^ ^7 — a in the g-plane. The intersection A\ fl A 2 — T is then the line 
5Re(g) = r = ^7 — a. 

The final transformation does, however, need tailoring to each domain separately as 
follows 

t = a(<7) := e ™ {2q -^ +2a)/2a q € Ai (A 11) 

£ = &( ? ) ;= e ^(29-7+2a)/( 7 r-2a) , £ ^ (A 12) 

so that the final composition transformations are 

£i(z) := exp^— [tan -1 w(z) - ^7 + a] J , 

£2(2) := expf j-^ - [tan" 1 w(z) - §7 + a] ) (A 14) 

where 



(A 13) 



g - cos 7 e' e - cos(2a - 7) 

W(z) = ; = — — -r . (A 15) 

sin 7 sin (7 — 2a) 

The image of A1/A2 is designed as the upper/lower half £-plane and V remains a shared 
boundary (see Figure 2). If we define f = £ + irj and $i(((x,y),r](x,y)) :— $i(x,y) 
(i = 1,2), the solutions for $1 and $2 are then available via Poisson's integral formula 
for the half plane 

•*■*>- i£^* 

The conditions (|A 4[) and (|A 5[) indicate that $1 and $2 are only non-zero on the image 
of r which is the positive real axis (i ^ 0) in the £-plane. The problem now boils down to 
determining the function f(z) :— u\ = uT 2 on T such that the stress matching condition 
(see (|2.7I) on V holds. Applying this condition is slightly non-trivial because the integrals 
(|A 16|) and (|A 17[) are formally singular for £ = ( + irj on T. They have well-defined 
(Cauchy principal) values by continuity with surrounding values of £ but taking normal 
derivatives of these integrals and subsequently computing them, nevertheless, requires 
due care. Consider the normal (77) derivative of $1 on V (t] = 0), for example. It is 
straightforward to show 



C-t 



(C-t) 2 + v 2 



dt. (A 18) 



and, after integration by parts, then 



Tl „ K ,o) - 1 f - 1 f ^Mtzp^m dt (A19) 
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since the Cauchy principal value of l/(i — Qdt is zero. The last integral on the right 
hand side of (|A 19j) is now regular. The symmetry of the velocity fields under y — ► — y 
in the z-plane can then be invoked to make the integration range finite. This refiectional 
symmetry carries over to the £-plane as the symmetry $j(l/t, 0) = $>i(t, 0) (i = 1,2) 
allowing, for example, l|A 16[) to be simplified to 



*i(C,»0 = - / *i(*,o) 

and (|XT9)) to 

$ M (C,o) = - / 



1 



(C-t) 2 +?? 2 (iC - I) 2 +* 2 /7 2 



dt. (A 20) 



d t+ iMm log ri-c 



$i, c (t,o)-$i, c (C,o) guM) 

t-C C* - 1 J ' TT '°\ ( 

(A21) 

These are the integral representations (along with the equivalent ones for <£> 2 ) used to 
impose the matching conditions and calculate the flow solution. 

In the matching process, the first step in determining / is to construct a global repre- 
sentation, f(9) = Yln=l c n^n(9), using 9 to parametrise T, c n as the expansion constants 
and the basis functions 

^n(S) '■= T2n{6 /8-max) ~ T2n-2{9 /Q-max)- (A 22) 

These are defined in terms of Chebyshev polynomials T n (9) := cos(n cos -1 9) with 
each designed to mirror the properties of /: f(±9 max ) = and df /d9\g—o — by the 
y— refiectional symmetry. This symmetry also means that the matching condition needs 
only to be applied (via collocation at the N positive zeros of T 2 n+i) over the upper 
half of r. It is tempting to carry out this procedure directly in the £— plane using the 
representation (|A 21j) and the sister integral for &2,rj- However, this proves inaccurate 
because both have an integrable singularity at t = (9 = ±9 max ). This causes loss of 
accuracy through two separate effects: a) the integrand has a singular derivative at t = 
so numerical quadrature is inefficient and b) the collocation points sparsely populate the 
neighbourhood of t = at extreme choices of a (— » or 7r/2) so the matching is not well 
imposed and convergence fails short of usual spectral (exponential) accuracy. Instead, 
the integral representations must be transformed to the physical z— plane and matching 
carried out there. 

The velocity profile along T is always smooth and typically only N — 20 or 30 is 
needed to see spectral drop off of 4-5 orders of magnitude. The limits a — * ir/2 and 
a — > 0, however, have to be treated carefully. For example, when a ^ 0.2 (» 10°) only 
a 100-pancl Simpson quadrature is needed to accurately calculate the integrals along T 
but this must be increased dramatically as a — > due to the extreme behaviour of the 
z = z(£) transformation in this limit (e.g. 10 4 panels proved sufficient for a — O(0. 001)). 
Once the solution is obtained, the fluxes Qi an d Q2 are calculated using Simpson's rule 
with typically 20 — 40 panels. This is the most costly part of the process as essentially a 
triple integral is being evaluated. Simple bisection in a is used to find a 'balanced' flux 
state where Q± + Q2 = for given /3, A and 7. 

As a final comment, it's worth remarking that the transformation q = q(z) (see the 
third subplot in figure [TJ]) achieves a separation of variables in the problem (the bound- 
aries are contours of constant r — Re(q)^. A solution could therefore be developed by 
separation of variables after a Fourier transform (in s) is taken of the inhomogeneity in 
the matching condition. The full procedure, however, boils down to essentially the same 



f The transformation q — q(z) is essentially a transformation to bipolar coordinates 
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problem of evaluating a triple integral albeit in this case the innermost one for u is an 
inverse Fourier transform and hence over a semi-infinite interval. 



Appendix B. Eccentric solutions 

The 'eccentric' solution has one fluid completely encapsulated by the other. For sake 
of argument, we describe the solution strategy for Ai in A\. The radius R and centre 
(er, 0) (with a < 0) define the geometry uniquely up to an arbitrary rotation around the 
duct axis and any reflection about a diameter neither of which, of course, affect the flux. 
The interface curve T is then 

r :={z\z = x + iy = a + Re lt> ; -tt < 6 tt } (B 1) 

which smoothly connects to the formula for T in the side-by-side solution (formally, R 
is +/— vc if r is convex/concave as viewed from x — -co: see the definition (|A 7[) ). The 
problem (|2.5[) — (I2.7[) is solved by conformally mapping the geometry of eccentric circles 
into one of concentric circles using a bilinear transformation £ = £(z). This is constructed 
by selecting a common pair of real inverse points (k, 0) and (is, 0) for T and the duct wall 
\z\ = 1 (so \kis\ = 1 and \k — a\\v — a\ = R 2 ) which ensures that the transformation 

Z — K , 

i-= B2 

z — v 

maps the two circles \z\ = 1 and T into concentric circles of radii (respectively) 

1-K R+(7-K 

za i := - & zv 2 ■= 7- (B3) 

1 — v R + a — is 

where 



K ) , = ±(l + f T 2 -i? 2 )-y/(l + ( 7 2 -fl 2 ) 2 -4^ 

(so v < — 1). In the £ = zue 1 ^ plane, the solution is found standardly using the expansions 



N I w 2n \ A + 1 

Hi = ^A n ^uj n - —^-J cosn<P + A Q log(ro/n7i) + — ^— (|z| 2 - 1), 



n=l 
N 



A — 1 , 



u 2 = B n w n cosncf) -\ — H- 2 ! 2 

n— 

which incorporate the boundary condition at \z\ = 1 (w = vj\) and the y— symmetry 
((f) — > —0) of the problem. The Fourier series in (/> of |z| = — «|/|^ — 1| and d\z\/dw 
on r need to be evaluated to apply the remaining matching conditions. This is done 
routinely using Simpson's rule with 200 panels when N — 100. In the limiting situations 
of (j — R. — > —1 (r approaching the duct wall) and a — > (approaching concentricity), 
these numbers are doubled to 400 and N = 200 to maintain at worst 10 -10 least square 
error in either matching condition. Calculation of the fluxes in A\ and A2 is again by 2D 
Simpson's rule using 100-200 panels per direction and simple bisection is used in R used 
to identify where Q\ + Q2 = for given (3, X and 1 + a — R (1 + a — R is fixed rather 
than a to avoid the complication of multiple solutions). 



REFERENCES 

Arakeri, J.H., Avila, F.E., Dada, J.M. & Tovar, R.O. 2000 Convection in a long vertical 



22 



R. R. Kerswell 



tube due to unstable stratification- A new type of turbulent flow? Current Science 79, 
859-866. 

Batchelor, G.K. & Nitsche, J.M. 1993 Instability of stratified fluid in a vertical cylinder. J. 

Fluid Mech. 252, 419-448. 
Beckett, F., Witham, F., Phillips, J.C. & Mader, H. private communication concerning 

experiments curently being carried out in the Department of Earth Sciences, University of 

Bristol - preprint coming. 
Charles, M.E. & Redberger, R.J. 1961 The reduction of pressure gradients in oil pipelines 

by the addition of water. Numerical analysis of stratified flows Can. J. Chem. Engng 40, 

70-75. 

Frigaard, LA. & Scherzer, O. 1998 Uniaxial exchange flows of Bingham fluids in a cylindrical 

duct IMA J. App. Math. 61, 237-266. 
Hasson, D., Mann, U. & Nir, A. 1970 Annular flow of two immiscible liquids. I. Mechanisms. 

Can. J. Chem. Engng. 48, 514. 
Huppert, H.E. & Hallworth, M.A. 2007 Bi-directional flows in constrained systems J. Fluid 

Mech. 578, 95-112. 

Joseph, D.D., Renardy, M. & Renardy, Y. 1984 Instability of the flow of two immiscible 
liquids with different viscosities in a pipe. J. Fluid Mech. 14, 309-317. 

Joseph, D.D., Nguyen, K. & Beavers, G.S. 1984 Non-uniqueness and stability of the con- 
figuration of flow of immiscible fluids with different viscosities. J. Fluid Mech. 14, 319-345. 

Joseph, D.D., Bai, R., Chen, K.P. & Renardy, Y.Y. 1997 Core-annular flows Ann. Rev. 
Fluid Mech. 29, 65-90. 

Lee, B.L. & White, J.L. 1974 An experimentalstidy of rheological properties of polymer 

melts in laminar shear flow and of interface defomration and its mechanisms in two-phase 

stratified flow. Trans. Soc. Rheol. 18, 467. 
Maclean, D.L. 1973 A theoretical analysis of bicomponent flow and the problem of interface 

shape Trans. Soc. Rheol. 17, 385. 
Markuskevich, A. I. 1965 Theory of Functions of a Complex Variable, vol 1 Prentice- Hall, 

Inc. Englewood Cliffs, New Jersey (p205-207) 
Minagawa, N. & White, J.L. 1975 Coextrusion of unfilled and TKh-filled polyethylene: influ- 
ence of viscosity and die cross-section on interface shape. Polymer Engng Sci. 15, 825. 
Moyers-Gonzalez, M.A. & Frigaard, LA. 2004 Numerical solution of duct flows of multiple 

visco-plastic fluids J. N on- Newtonian Fluid Mech. 122, 227-241. 
Seon, T, Znaien, J., Salin, D., Hulin, J. P., Hinch, E.J. & Perrin, B. 2007 Transient 

buoyany-driven front dynamics in nearly horizontal tubes Phys. Fluids 19, 123603 
Southern, J.H. & Ballman, R.L. 1973 Stratified bicomponent flow of polymer melts in a 

tube Appl. Polymer Symp. 20, 175-189. 
Taghavi, S.M., Seon, T., Martinez, D.M. & Frigaard, LA. 2009 Buoyancy-dominated 

displacement flows in near-horizontal channels: the viscous limit J. Fluid Mech. 639, 1-35. 
Vlasov, V.I. 1986 Solution of a Dirichlet problem in a crescent-shaped domain J. Eng. Phys. 

& Thermophys. 50, 741-747. 
White, F. M. Viscous Fluid Flow McGraw-Hill (pl24) 

Williams, M.C. 1975 Migration of two liquid phases in capillary extrusion: an energy interpre- 
tation AIChE. J. 21, 1204. 

Yu, H.S. & Sparrow, E.M. 1967 Strained laminar flow in ducts of arbitrary shape AIChE. J. 
13, 10. 

Znaien, J., Hallez, Y., Moisy, F., Magnaudet, J., Hullin, J. P., Salin, D. & Hinch, E.J. 
2009 Experimental and numerical investigations of flow structure and momentum transport 
in a turbulent buoyancy-driven flow inside a tilted tube. Phys. Fluids 21, 115102. 



